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ABSTRACT 

We investigate a novel Bayesian analysis method, based on the Stochastically Lighting Up Galaxies 
(slug) code, to derive the masses, ages, and extinctions of star clusters from integrated light pho¬ 
tometry. Unlike many analysis methods, slug correctly accounts for incomplete IME sampling, and 
returns full posterior probability distributions rather than simply probability maxima. We apply our 
technique to 621 visually-confirmed clusters in two nearby galaxies, NGG628 and NGC7793, that 
are part of the Legacy Extragalactic UV Survey (LEGUS). LEGUS provides Hubble Space Telescope 
photometry in the NUV, U, B, V, and I bands. We analyze the sensitivity of the derived cluster 
properties to choices of prior probability distribution, evolutionary tracks, IME, metallicity, treat¬ 
ment of nebular emission, and extinction curve. We find that slug’s results for individual clusters are 
insensitive to most of these choices, but that the posterior probability distributions we derive are often 
quite broad, and sometimes multi-peaked and quite sensitive to the choice of priors. In contrast, the 
properties of the cluster population as a whole are relatively robust against all of these choices. We 
also compare our results from slug to those derived with a conventional non-stochastic fitting code, 
Yggdrasil. We show that slug’s stochastic models are generally a better fit to the observations than 
the deterministic ones used by Yggdrasil. However, the overall properties of the cluster populations 
recovered by both codes are qualitatively similar. 

Subject headings: methods: data analysis — methods: statistical — galaxies: individual (NGG628, 

NGG 7793) — galaxies: star clusters: general — techniques: photometric 


1. INTRODUCTION 

Star clusters represent a scale of star formation in¬ 
termediate between individual stars and entire galaxies. 
They likely represent the gravitationally-bound peaks of 
a continuous distribution of star formation across length 
and time scales. Eor this reason, the study of their prop- 
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erties, and any variation in those properties as a function 
of galactic environment, is critical to an overall theo¬ 
retical understanding of the star formation process (see 
Portegies Zwart et al. 2010, Longmore et al. 2014, and 
Krumholz 2014 for recent reviews). 

While the study of star clusters is an old topic, except 
in the Milky Way and the very closest external galax¬ 
ies, much of this work has by necessity focused on the 
relatively massive clusters, ^ 10^’^ Mq or more (e.g. 
Zhang & Eall 1999; Larsen & Richtler 2000; Bik et al. 
2003; de Grijs et al. 2003b; Eall et al. 2005; Goodwin & 
Bastian 2006; Larsen 2009). This has created significant 
challenges for theoretical interpretation, for two reasons. 
The first is simply statistics: as discussed in more detail 
below, observed star cluster mass functions appear to be 
well-fit by powerlaws of the form dN/dM oc M~^ (e.g., 
Williams & McKee 1997; Zhang et al. 1999; Larsen 2002; 
Bik et al. 2003; de Grijs et al. 2003b; Goddard et al. 
2010; Bastian et al. 2011; Eall & Ghandar 2012; Eoues- 
neau et al. 2012). Such a mass function implies that clus¬ 
ters in the mass range 10^’^ — 10^’^ Mq are nearly ten 
times as numerous as those in the range ^ 10^’^ — 10^’^ 
Mq that is typically studied, so a sample consisting only 
of massive clusters necessarily has much less statistical 
power than one probing to lower masses. The second 
is that theoretical models for the evolution and disrup¬ 
tion of star clusters, either due to stellar feedback or due 
to environmental influence, make some of their strongest 
predictions for the effects of cluster mass at masses be¬ 
low this range (e.g. Earners et al. 2005; Parmentier et al. 
2008; Eall et al. 2010; Kruijssen 2012; Kruijssen et al. 
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2012 ). 

In the last five years a number of groups have made 
a concerted effort to extend the star cluster data set to 
lower masses and a broader range of galactic environ¬ 
ments, a necessary step if we are to differentiate between 
models. Doing so is challenging both observationally 
and theoretically. Observationally, ground-based stud¬ 
ies probing to lower masses are mostly limited by reso¬ 
lution to the Milky Way (Borissova et al. 2011) and the 
Magellanic Clouds (Hunter et al. 2003; Rafelski & Zarit- 
sky 2005). Extending this sample is one of the primary 
goals of two ongoing Hubble Space Telescope programs: 
the Panchromatic Hubble Andromeda Treasury (PHAT; 
Dalcanton et al. 2012) and the Legacy Extra-Galactic 
UV Survey (LEGUS; Galzetti et al. 2015). The former 
focuses on M31 to extreme depth, providing an extensive 
and very complete cluster catalog, but for only a single 
galaxy. The latter targets 50 nearby galaxies chosen from 
across the star-formiong portion of the Hubble sequence, 
providing data on a much greater range of galactic en¬ 
vironments, but with less depth. A preliminary cluster 
catalog for PHAT is now available (Johnson et al. 2012), 
and catalogs for the LEGUS galaxies will be released in 
the future (Adamo et ah, 2015, in preparation). 

The theoretical challenge when working with these ex¬ 
tended cluster samples is that deriving masses and other 
properties for small clusters is non-trivial because they 
do not fully sample the stellar initial mass function (e.g., 
Gerviho & Vails-Gabaud 2003; Cerviho & Luridiana 
2004, 2006; Mafz Apellaniz 2009). Traditional methods 
of determining star cluster properties by isochrone fit¬ 
ting therefore begin to fail, because at such small masses 
the relationship between clusters’ photometric properties 
and their mass, age, and other physical properties is non- 
deterministic. One approach to dealing with this prob¬ 
lem is to resolve at least the most massive members of the 
cluster and determine their properties star-by-star using 
a color-magnitude diagram (CMD, Beerman et al. 2012). 
However, for young open clusters where the stellar sur¬ 
face density is high and confusion is significant, this re¬ 
quires extreme resolution and depth in the observations, 
and thus is not yet practical as a method for obtain¬ 
ing large cluster samples across a wide range of galactic 
environments.^^ The alternative approach is to develop 
new statistical techniques that can cope with partial sam¬ 
pling, and this has motivated the development of three 
major stochastic stellar population synthesis and analysis 
codes: MASSclean (Popescu & Hanson 2009, 2010a, b), a 
stochastic version of pegase (Eouesneau & Lancon 2010; 
Eouesneau et al. 2012, 2014), and slug (da Silva et al. 
2012, 2014; Krumholz et al. 2015). While these codes use 
slightly different statistical techniques, they all attempt 
to solve essentially the same problem: given a set of ob¬ 
served photometric properties for a star cluster, what 
should we infer about the probability distribution for its 
mass, age, extinction, or other physical properties? 

In this paper we use slug, the Stochastically Lighting 

To our knowledge the largest-scale application of this tech¬ 
nique to appear in the literature to date is in Fouesneau et al. 
(2014), who use PHAT photometry to obtain CMD-based ages for 
iOO clusters in M51. However, this paper presents the CMD-based 
results only for the purposes of comparison with those based on 
integrated photometry, and does not contain the CMD analysis 
itself, which is deferred to an as-yet unpublished paper. 


Up Galaxies code, and its post-processing tool for anal¬ 
ysis of star cluster properties, cluster_slug, to analyze 
an initial sample of clusters from the Legacy Extragalac- 
tic Galaxy Ultraviolet Survey (LEGUS; Galzetti et al. 
2015). The clusters were identified in three HST WEC3 
pointings (field of view 2.'7 x 2.'7). Two pointings span 
the radius of the inner disk of NGC7793, and include 
its nucleus, while the third pointing is in NGC628, just 
inside of the eastern inner disk edge. Both NGG 628 and 
NGC 7793 are well-studied late-type spiral galaxies, and 
are relatively isolated relative to other massive systems. 
However, they provide complementary views of star clus¬ 
ter populations. Whereas NGC 7793 is in the nearby 
Sculptor group at 3.6 Mpc, and is inclined, NGC 628 is 
nearly face-on and is more distant, at 9 Mpc (Tully et al. 
2009). Their star formation rates and stellar masses also 
differ by a factor of ^ 3, and are 0.7 Mq yr“^ and 4 x 10^ 
Mq for NGC 7793, and 2.0 M© yr-^ and 1.6 x 10^^ M© 
for NGC628 (Lee et al. 2009; Cook et al. 2014). Thus 
NGC 628 and 7793 represent examples towards the ex¬ 
tremes of the full LEGUS sample of late-type spirals in 
terms of distance, orientation, mass, and star formation 
rate. This makes them a useful testbed for the perfor¬ 
mance of our analysis, one that will be extended to the 
full LEGUS sample in future work. 

We have two goals in this paper. Eirst, we wish to un¬ 
derstand the performance of the stochastic stellar popu¬ 
lation analysis, including how it compares to traditional 
fitting methods, and to understand any systematic er¬ 
rors that might appear in the results due to choices of 
stellar tracks, stellar IME, extinction law, assumed pri¬ 
ors, or similar issues. Analysis of this sort has previously 
been performed by Popescu et al. (2012) for the Large 
Magellanic Cloud, by Fouesneau et al. (2012, 2014) for 
M31, and by de Meulenaer et al. (2013, 2014, 2015) for 
M31 and M33, but studies of systematics have mostly 
been limited to questions of isolated metallicity, evolu¬ 
tionary tracks, and extinction law. To our knowledge no 
previous work has covered the possible systematics thor¬ 
oughly in the context of a single framework, and none 
at all have explored the issue of priors. Second, we wish 
to describe a portion of the LEGUS data pipeline, which 
will be used to produce star cluster catalogs based on 
cluster_slug modeling for the entire LEGUS sample. 
While this paper focuses on three fields in two galax¬ 
ies for early analysis, the methods we develop here will 
be applied to the entire LEGUS sample. Ultimately, we 
will provide a catalog of homogeneously-observed and - 
analyzed star clusters, whose properties have been de¬ 
rived with a fully stochastic treatment of stellar popu¬ 
lation synthesis, for 50 galaxies that sample across the 
star-forming part of the Hubble sequence. This data set 
will provide tremendous insight into the extent to which 
star cluster formation and evolution depends on galactic 
environment. 

The plan for the remainder of this paper is as follows. 
In Section 2 we describe the method we use to derive the 
physical properties of clusters, both stochastically and, 
for comparison, using a deterministic method. Section 3 
presents the results of this analysis for individual clus¬ 
ters, and discusses general trends and possible systemat¬ 
ics. In this section we also investigate how the stochastic 
and deterministic analyses compare. Section 4 extends 
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this analysis to the properties of the star cluster pop¬ 
ulation as a whole. Finally, Section 5 summarizes our 
conclusions. 

2. METHODOLOGY 

2.1. The Input Photometric Catalog 

A description of the steps required to produce final 
cluster catalogues of the LEGUS targets can be found in 
Galzetti et al. (2015), and in Adamo et al. (in prep) we 
will present the custom pipelines we have developed to 
analyze the whole sample in a homogeneous fashion. We 
summarize here the main aspects of this process which 
led to the final cluster catalogues of the three pointings in 
two galaxies used in this work, i.e., NGC628 East (here¬ 
after, NGG628e), NGG 7793 West (NGC7793w), and 
NGG 7793 East (NGG 7793e). We create our catalogs by 
using the SExtractor algorithm (Bert in & Arnouts 1996) 
on white-light images produced with the five standard 
HST LEGUS UV - UBVI bands (WFG3 F275W and 
F336W, WFG3 or ACS F438W or F435W, F555W, and 
F814W, respectively). The exact filters and exposure 
times used in each pointing are listed in Table 1. The 
cluster candidate catalogues contain only sources which 
satisfy the following criteria: 1) the V band concentra¬ 
tion index (Cl) must be greater than the stellar Cl peak 
(the Cl reference value slightly changes as function of 
galactic distance and HST camera); 2) the cluster can¬ 
didate must be detected in two contiguous bands with a 
signal-to-noise higher than 3. 

The multi-band photometry of the cluster candidates 
of the three galactic fields used in this work has been 
derived with a fixed photometric aperture. We used an 
aperture radius of 4 and 5 px for NGC 628 and NGC 7793, 
respectively. The aperture radius is determined to in¬ 
clude at least 50% of the flux of a median built growth 
curve of several isolated clusters selected in the field of 
each galaxy. The value will change mainly as function 
of galactic distance. The sky annulus 1 pixel wide is lo¬ 
cated at 7 pixel from the centre of the source. Averaged 
aperture corrections in all the filters have been estimated 
using manually selected isolated clusters in each galac¬ 
tic field. To take into account uncertainties produced 
by using a fixed aperture correction, we have added in 
quadrature the standard deviations of the aperture cor¬ 
rection and the raw photometric error. 

To remove contaminants from these automatically cre¬ 
ated catalogues we have visually inspected a subsample 
of cluster candidates in each galactic field with detections 
in at least 4 filters and absolute magnitudes brighter 
than My = —6 mag. Each inspected source has then 
been assigned a morphology quality flag. Class 1 objects 
are compact and centrally concentrated clusters; class 2 
systems are clusters which show slightly elongated den¬ 
sity profiles; class 3 are systems which show asymmetric 
profiles and multiple peaks, suggesting an association of 
stars; class 4 objects are single stars or artifacts, or any 
other spurious detection which can be excluded from the 
cluster catalogue. In this work we will only use high- 
confidence, visually-confirmed cluster catalogues clusters 
with photometry in at least 4 filters and flagged as class 

1,2, and 3. However, as a service for readers who might 
wish to make their own classifications and produce their 
own samples, in the electronic edition of the paper we 


provide a separate machine-readable table containing the 
results of a cluster_slug analysis of the clusters we omit 
here. We also note that, given their morphology, one 
might question whether class 3 objects should be called 
“clusters”. For the purposes of this analysis we treat 
them as such, but the fact that we have included them 
should be kept in mind when examining our results on 
the cluster population. Excluding them would produce a 
population that includes somewhat fewer young objects, 
since we find slightly younger than average ages for class 
3 objects. However, since the focus of this work is on the 
method of analysis rather than the details of the popu¬ 
lation, we defer further discussion of this topic to future 
work. 

Before proceeding, we make two additional notes re¬ 
garding the catalog. First, we exclude the nuclear star 
cluster of NGG 7793 from all the analysis presented be¬ 
low. Although this cluster is flagged to be in the high- 
confidence catalog, it is almost certainly not a simple 
stellar population, and we should not analyze its prop¬ 
erties as if it were. Second, due to partial overlap of the 
NGG 7793e and NGG 7793w pointings, some clusters ap¬ 
pear twice in the catalog, once for each field. We have 
retained the two separate entries in the analysis below, 
since they are observed with slightly different filters and 
thus do not produce completely identical results. How¬ 
ever, when presenting population statistics in Section 4, 
we include only the version of the cluster that appears in 
our NGG 7793e catalog. 

The final, high-confidence cluster catalog we use for 
this paper contains 621 distinct clusters. Gounting 
the duplicates that appear in the overlapping fields for 
NGG 7793 separately, the total rises to 645. 

2.2. Stochastic Models with cluster_slug 

Our stochastic analysis makes use of the 
cluster_slug software package that is part of the 
slug (Stochastically Lighting Up Galaxies) software 
suite (da Silva et al. 2012, 2014; Krumholz et al. 2015). 
The Bayesian analysis performed by cluster_slug takes 
as input a set of absolute magnitudes Mp in some set 
of filters, with corresponding errors AMp. It returns 
as output the posterior probability distribution for the 
cluster mass M, age T, and extinction Ay. Details 
of how this operation is performed can be found in 
the papers referenced above, and both the slug suite 
and the full software pipeline used in this paper are 
available at https : //www. slugsps . com. Here we limit 
our discussion to the choice of two key inputs for the 
analysis: the libraries of simulated star clusters on 
which cluster_slug operates and the choice of prior 
probabilities that enter any Bayesian analysis. For a 
third input, the kernel density estimation bandwidth, 
we use h = 0.1 dex in the physical variables and h = 0.1 
mag in the photometric ones. 

The cluster_slug code requires a library of simulated 
star clusters to act as a “training set”. We produce these 
using the slug stochastic stellar population synthesis 
code. To generate a library, we must specify the stel¬ 
lar evolution tracks, metallicity, extinction curve, stel¬ 
lar initial mass function (IMF), and the fraction of the 
ionizing light that is reprocessed into nebular emission 
within the observational aperture. One of the goals of 
this paper is to study how these choices affect the de- 
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Table 1 

List of filters and exposure times 


Field Filter tgxp (s) Filter texp (s) Filter tgxp (s) Filter texp (s) Filter tgxp (s) 


NGC628e WFC3 F275W 2361 WFC3 F336W 1119 ACS F435W 4720 WFC3 F555W 965 ACS 814W 1560 

NGC7793W WFC3 F275W 2349 WFC3 F336W 1101 WFC3 F438W 947 WFC3 F555W 680 WFC3 F814W 430 

NGC7793e WFC3 F275W 2349 WFC3 F336W 1101 WFC3 F438W 947 ACS F555W 1125 ACS F814W 971 



log T [yr] 


Figure 1. Absolute Vega magnitude versus cluster age for Padova 
tracks with Z = 0.02 and 0 = 0.5, Geneva tracks with Z = 0.014 
and (j) = 0.5, and Padova tracks with Z = 0.02 and 0 = 0, in 
the filters WFC3 F275W, WFC3 F555W, and WFC3 F814W, 
as predicted by slug for a 10® Mq cluster with a fully-sampled 
Kroupa IMF. 

duced cluster properties, and for this reason we generate 
a range of libraries, summarized in Table 2. The param¬ 
eters we vary in these libraries are as follows: for tracks 
we consider both those from the Padova group including 
TP-AGB stars (Vassiliadis & Wood 1993; Girardi et al. 
2000; Vazquez & Leitherer 2005) and from the Geneva 
group (Ekstrom et al. 2012). For metallicities we con¬ 
sider Z = 0.004, 0.008, 0.020 for the Padova tracks, and 
Z = 0.014 for the Geneva tracks.We consider extinc¬ 
tion curves appropriate to both the Milky Way (Landini 
et al. 1984; Fitzpatrick 1999) and to a starburst (Calzetti 
et al. 2000), and both Kroupa (2001) and Ghabrier (2005) 
IMFs. Finally, for nebular emission we consider both 
0 = 0.5 (50% of ionizing photons available to produce 
nebular emission inside the aperture) and 0 = 0 (no neb¬ 
ular emission). Figure 1 provides examples of how some 
of these choices affect the predicted variation of magni¬ 
tude versus age in select filters. All the examples shown 
are for a 10^ Mq cluster with a fully-sampled Kroupa 
IMF, and no extinction. One can clearly see in the Fig¬ 
ure 1 the differences in predicted luminosity due to the 
treatment of TP-AGB stars at older ages, and due to the 
treatment of nebular emission at younger ages. 

In principle it would be best to treat metallicity as a continu¬ 
ous variable to be fit, as we will fit cluster mass and age. However, 
tracks including TP-AGB stars that are sufficiently finely sampled 
in metallicity to make such a treatment possible did not become 
available from the Padova group until late 2014, and have not yet 
been incorporated into slug. No correspondingly fine tracks are 
available from the Geneva group yet. As a result, we will treat 
metallicity as fixed. As discussed below, this is not a significant 
limitation for our target galaxies, because the expected metallicity 
variation within the galaxy is relatively small. 


Our fiducial library, which we use for all purposes 
unless specified otherwise, is pad_020_kroupa_MW. This 
uses Padova isochrones, metallicity Z = 0.02, a Kroupa 
IMF, a Milky Way extinction curve, and 50% of ionizing 
photons converted to nebular emission. We adopt this 
as our fiducial model because it appears to be the closest 
match to what we know of the target galaxies. Observed 
IMFs in local spiral galaxies are consistent with a univer¬ 
sal function consistent with the Kroupa (2001) parame¬ 
terization (see the recent reviews by Krumholz (2014) 
and Offner et al. (2014)). Both NGC628 and NGC 7793 
have mean metallicities close to Solar, and for the typ¬ 
ical ^ 0.03 dex kpc“^ metallicity gradients observed in 
the inner disks of local spirals (e.g., Pilyugin et al. 2004), 
the < 10 kpc-wide regions spanned by our observed fields 
should have end-to-end metallicity differences of only a 
few tenths of a dex. On average we observe that, for 
the deterministic models at least (see below), clusters 
within face-on Solar metallicity spirals are better fitted 
by a Milky Way extinction curve (Cardelli et al. 1989) 
than with the alternatives. Finally, we use the Padova 
rather than Geneva evolutionary tracks because the ma¬ 
jority of our clusters are older than ^ 10 Myr, and in 
this age range the superior treatment of AGB stars in 
the Padova models is advantageous. All these fiducial 
choices are discussed further in Adamo et al. (2015, in 
preparation). 

In addition to these choices, for each library we must 
specify how the masses, ages, and extinctions of the sim¬ 
ple stellar populations are to be sampled. For all the li¬ 
braries used in this paper, we choose a distribution that 
varies with cluster mass and age as 

P\ih{x) = PMilogM)pT{\ogT)pAv{Av) (1) 


where 




PM(logM)oc < 

fr 

1 20“dogM-4)^ 

2 < log M < 4 

4 < log M < 8 

(2) 

PT(iogr)oc < 

fr 

Q0-(log'r-8), 

5<logT<8 

8 < log T < log Tmax 

(3) 


PAv-ocl, 0 < Av-< 3 (4) 

In the above expressions M is in units of Mq, T is in 
units of yr, and Ay is in mag. We use T^ax = 15 Gyr for 
models using the Padova tracks, and T^ax = 1 Gyr for 
models using the Geneva tracks. This sampling distribu¬ 
tion places most realizations at small masses and ages, 
where stochastic variation is largest, and uses fewer real¬ 
izations at higher masses and older ages where they are 
smaller. We generate 10^ realizations for each library. 

The final ingredient needed for a cluster_slug calcu¬ 
lation is a prior probability distribution for the physical 
properties, which we use to weight the simulations in the 
library. One can demonstrate that the choice of priors 
will be important in at least some regimes via a simple 
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Table 2 

Model libraries for cluster.slug 


Name 

Tracks 

IMF 

Z 

Extinction^ 


log M 
(Me) 

logT 

(yr) 

(mag) 

^ Realizations 

pad_020_kroupa_MW^ 

Padova 

Kroupa 

0.020 

MW 

0.5 

2 

- 8 

5 - 10.18 

0-3 

10'^ 

pad_020_kroupa_SB 

Padova 

Kroupa 

0.020 

SB 

0.5 

2 

- 8 

5 - 10.18 

0-3 

10'^ 

pad_004_kroupa_MW 

Padova 

Kroupa 

0.004 

MW 

0.5 

2 

- 8 

5 - 10.18 

0-3 

10'^ 

pad_008_kroupa_MW 

Padova 

Kroupa 

0.008 

MW 

0.5 

2 

- 8 

5 - 10.18 

0-3 

10^ 

pad_020_chabrier_MW 

Padova 

Chabrier 

0.020 

MW 

0.5 

2 

- 8 

5 - 10.18 

0-3 

10'^ 

gen_014_kroupa_MW 

Geneva 

Kroupa 

0.014 

MW 

0.5 

2 

- 8 

5- 9.00 

0-3 

10'^ 

pad_020_kroupa_MW_noneb 

Padova 

Kroupa 

0.020 

MW 

0.0 

2 

- 8 

5 - 10.18 

0-3 

10'^ 


^ MW = Milky Way extinction curve, SB = starburst extinction curve 

^ ^ is the fraction of the ionizing photons that produce nebular emission within the aperture; it combines the effects of a 
covering fraction < 1 and some portion of the ionizing photons being absorbed directly by dust 
^ Fiducial model 


thought experiment. Consider what might seem like a 
natural prior, a distribution that is flat in the log of the 
age. Stars’ luminosities and colors do not change sig¬ 
nificantly until their ages reach ^ 1 Myr. Thus all star 
clusters younger than this age look identical, and the 
shape of the posterior probability distribution at ages 
below ^ 1 Myr must therefore be the same as the shape 
of the prior probability distribution. For a prior that 
is flat in log age, we would therefore conclude that age 
ranges of log(T/yr) = 3.5 — 4 and log(T/yr) = 5.5 — 6 are 
equally likely. This seems unlikely to be correct: even if 
it were completely unbound, a stellar population formed 
< 1 Myr ago would still be close enough together to be 
recorded as a cluster in the LEGUS catalog, so every 
cluster that reaches an age of log(T/yr) = 3.5 — 4 must 
eventually reach an age of log(T/yr) = 5.5 — 6. However, 
this cluster will reside in the interval log(T/yr) = 3.5 — 4 
for a mere 6,800 yr, while it will spend 680,000 yr in 
the interval log(T/yr) = 5.5 — 6. Clearly there should 
be 100 times as many clusters in the latter bin as in 
the former, making the latter 100 times as likely. Even 
worse, suppose further that we have a cluster whose col¬ 
ors admit ages older as well as younger than ^ 1 Myr. 
In this case the weight that we would end up assigning 
to the larger ages would depend strongly on whether our 
model grid contained a youngest age of, say, 10^ yr, 10^ 
yr, 10^ yr, or 10^ yr, because the amount of probability 
“phase space” allowed at young ages would depend on 
this choice. Clearly some attention to priors is required 
to avoid outcomes such as this. 

For ages younger than a few Myr, the proper prior is 
almost certainly fiat in age. This is because, even if a 
newly-formed collection of stars has negligible gravita¬ 
tional binding energy, at the typical velocity dispersions 
of a few km s“^ found in young clusters, it will disperse 
over a region no more than ^ 10 pc in size over this time. 
In an extragalactic survey such as LEGUS, it would still 
be detected as a cluster. Thus cluster dispersal is ir¬ 
relevant at such young ages. Similarly, as noted above, 
stellar evolution effects are also negligible. Given the 
absence of any physical mechanism to bias the age dis¬ 
tribution, all ages are equally likely, suggesting that the 
proper prior is fiat in age. 

For ages greater than ^ 10^’^ yr, the proper prior is sig¬ 
nificantly less certain, as there is an extensive debate in 
the literature on this topic. Some authors report distri¬ 
butions dNIdT oc T~^ (i.e., flat in log age; e.g. Fall et al. 


2005, 2009; Chandar et al. 2010a, b, 2011; Fall & Chan- 
dar 2012; Fouesneau et al. 2012; Popescu et al. 2012), 
while others (mostly focusing on ages > 10 Myr) find 
dN /dT oc (i.e., fiat in age; e.g. Boutloukos & Earners 
2003; de Grijs et al. 2003a; Gieles et al. 2007; Bastian 
et al. 2011, 2012b, a; Silva-Villa et al. 2014); some authors 
also report intermediate results, with the index of the age 
distribution changing as a function of age (e.g. Fouesneau 
et al. 2014). The results appear to depend in part on the 
criteria one uses to define what is a cluster and thus 
should be included in the cluster catalog; see Krumholz 
(2014) for a recent review. The situation is further com¬ 
plicated by the fact that the prior we want to use is not 
the intrinsic distribution of star cluster ages, but the con¬ 
volution of that with our selection function. Properly 
modeling this would require much greater understand¬ 
ing of our observational completeness and biases than 
we now possess. Given these complexities, as a fiducial 
prior we choose a compromise, dN/dT ^ which is 

intermediate between flat in age {dN/dT ^ constant) 
and flat in log age {dN/d log T ^ constant, so that 
dN/dT ^ T~^). However, we will explore the effects 
of varying this choice below. The maximum possible age 
will be set by the maximum age in our library. 

For the prior on the mass distribution, we note that 
observations of young star clusters consistently find 
dN/dM oc (e.g., Williams & McKee 1997; Zhang 
et al. 1999; Larsen 2002; Bik et al. 2003; de Grijs et al. 
2003b; Goddard et al. 2010; Bastian et al. 2011; Fall & 
Ghandar 2012; Fouesneau et al. 2012), with relatively lit¬ 
tle variation. These results are mostly derived at masses 
high enough that stochasticity is unimportant, and thus 
it seems reasonable simply to extrapolate them into our 
regime. We therefore adopt p(logM) oc M~^ (corre¬ 
sponding to dN/dM oc M“^) as our fiducial prior on the 
mass distribution. Two caveats are worth mentioning, 
though. First, as with the age distribution, the truly cor¬ 
rect prior to use would be the convolution of this function 
with our observational selection, which unfortunately is 
not fully characterized as yet. Second, to avoid a log¬ 
arithmic divergence, we must truncate the prior distri¬ 
bution at both the high and low mass ends; we implic¬ 
itly choose the truncation values via the range of masses 
we include in our library. Fortunately the choices here 
have relatively little effect, because our library covers 
masses from 10^ — 10^ Mq, which spans the entire plan- 
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sible mass range: the lower limit of 100 Mq is almost 
certainly smaller than the smallest cluster we have any 
hope of detecting, and, given the size of our sample, the 
upper limit is much more massive than the most massive 
cluster we would expect to find even if the underlying 
cluster mass function continued as dNjdM oc M~^ to 
even higher masses. 

Finally, for the Ay distribution, observations indicate 
that Ay values tend to be at most ^ 0.5 at ages above 
^ 10 Myr (e.g., Whitmore & Zhang 2002; Whitmore 
et al. 2011; Mengel et al. 2005; Bastian et al. 2014; Hol- 
lyhead et al. 2015). At younger ages, Fouesneau et al. 
(2012, 2014) report a fairly broad distribution of Ay. In 
principle we could use a prior that combines age and Ay 
to disfavor a combination of old age and high-Ay. How¬ 
ever, the functional form of this constraint is not well- 
characterized, and we wish to keep our priors on age and 
Ay independent to ease the analysis of how our results 
depend on a choice of priors. For this reason we choose 
to adopt a flat distribution in Ay as a prior. 

In summary, our fiducial prior distribution follows 
Pprior{x) OC with no dependence on Ay. How¬ 

ever, in Section 3.3 we will check the dependence of our 
results upon this choice. 

2.3. Deterministic Models and Traditional -Fitting 

Method 

The standard photometric analysis of the LEGUS clus¬ 
ter catalogues is performed with Yggdrasil^^ (Zack- 
risson et al. 2011) deterministic models and a tradi¬ 
tional x^-fitting procedure (e.g., see Calzetti et al. 2015). 
Yggdrasil is a spectral synthesis code that can provide 
spectral and integrated flux tables for a large set of stel¬ 
lar and nebular parameters. In Yggdrasil, the spectrum 
produced by the stellar population at each age step is 
then used as input for a calculation of the nebular emis¬ 
sion using cloudy (Ferland et al. 2013) to ensure a self- 
consistent treatment. 

For the LEGUS cluster analysis in this paper, we 
have carefully matched the physical models used in 
Yggdrasil to those of the fiducial cluster_slug library. 
Thus for our Yggdrasil models we assume that clusters 
are a simple stellar population, and we compute the stel¬ 
lar spectrum using Starburst99 (Leitherer et al. 1999; 
Vazquez & Leitherer 2005) atmosphere models. We use 
the Kroupa (2001) IME, the same Padova AGB stellar 
tracks used in the cluster_slug models, and a metallic- 
ity Z = 0.02. We calculate nebular emission assuming 
that (1) the metallicity of gas and stars is the same, (2) 
the average gas density is 100 cm“^, (3) the ionizing flux 
is normalized by setting it to that of a 10^ Mq stel¬ 
lar population, and (4) the covering factor set to 50%, 
meaning that 50% of the ionizing photons are assumed 
to be absorbed by hydrogen atoms within the observa¬ 
tional aperture. To model extinction, we adopt a Milky 
Way extinction curve, and consider differential extinc¬ 
tions E(B—V) between 0 and 1.5 mag with steps of 0.01 
mag. We apply extinction to the model spectra before 
convolving them with the filter throughput. Thus in ev¬ 
ery respect except nebular emission, which we discuss 

http://ttt.astro.su.se/proj ects/yggdrasil/yggdrasil. 

html 


in Section 2.3.1, our Yggdrasil models are completely 
matched to our fiducial cluster_slug one. 

The x^-fitting algorithm used to derive cluster physical 
properties has been presented and tested in Adamo et al. 
(2010). The algorithm compares the observed photome¬ 
try to the library of Yggdrasil models, and returns the 
library cluster age, extinction, and mass that deviates 
least from the observations. A quality of the fit value is 
estimated from the reduced x^, he., Q-parameter. Good 
fits have Q values close to 1, while values below 0.1 sug¬ 
gest poor fits to the observed cluster SEDs. The x^- 
fitting algorithm has also been extended by Adamo et al. 
(2012) to produce uncertainties in the derived cluster 
physical properties contained within the 68% confidence 
limits of the solution parameter space (Galzetti et al. 
2015, their Eigure 9). 

2.3.1. A Caveat Regarding Nebular Emission 

Before presenting results, we pause to add a caveat 
to our comparison of the results obtained with 
Yggdrasil and duster_slug. We have tried to ensure 
that the underlying physical models are as similar as pos¬ 
sible in the two codes, so that any differences between 
them arise purely from the fact that cluster_slug in¬ 
cludes stochasticity and returns a full posterior PDE, 
while our models using Yggdrasil are non-stochastic, 
and our fits to them use a traditional x^ method. The 
one quantity we have not been able to match precisely is 
the treatment of nebular emission. 

To good approximation, an H ll region produces a neb¬ 
ular luminosity per ionizing photon injected that is a 
function of the gas metallicity and ionization parame¬ 
ter. While the Yggdrasil and cluster_slug models are 
matched in metallicity, they are not precisely matched in 
ionization parameter. As noted above, the nebular con¬ 
tribution in the Yggdrasil models we use in this paper 
has been derived by taking spectra computed for a 10^ 
Mq stellar population illuminating an H ll region with 
a constant density of 10^ cm“^ and a 50% covering frac¬ 
tion. (See Zackrisson et al. (2011) for further details.) 
In the Yggdrasil models, the inner radius of the H II 
region scales with total luminosity to the 1/2 power, so 
the ionization parameter at the illuminated face of the 
H II region is kept roughly constant - roughly rather than 
precisely because the ratio of ionizing to total luminosity 
is not constant as a stellar population ages. However, 
because the volume of ionized gas changes as the lumi¬ 
nosity does, and the ionizing photon flux changes as the 
photons propagate through the nebula, this prescription 
causes the volume-averaged ionization parameter to vary 
significantly with the mass and age of the stellar popu¬ 
lation. The mean is close to log U = —2.5, but there is 
a non-negligible spread. 

Slug’s treatment of nebular emission also assumes con¬ 
stant density within an H II region, and we have chosen 
the same 50% covering fraction as in the Yggdrasil mod¬ 
els. However, in its simplest form, slug assumes that 
the volume-averaged ionization parameter has a con¬ 
stant value log U = —3, while the ionization param¬ 
eter at the inner edge is not held constant.^^ While 

Full details of the slug method, and the computational moti¬ 
vations for choosing a constant ionization parameter, are given in 
Krumholz et al. (2015). 
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slug also supports an option to pass the spectra through 
cloudy, and thus in principle we could fully emulate 
Yggdrasil’s assumptions, it is not computationally fea¬ 
sible to run cloudy on all, or even a substantial sub¬ 
set, of the nearly 10^ models in our libraries. Thus we 
are left with an imperfect match in ionization parame¬ 
ter between the cluster_slug and Yggdrasil models; 
the Yggdrasil models assume that H II regions form a 
sequence where the ionization parameter at the illumi¬ 
nated face is fixed but the volume-averaged ionization 
parameter is not, while the duster_slug models fix the 
volume-averaged but not the illuminated face ionization 
parameters. 

How does this difference in ionization parameter trans¬ 
late into a difference in photometry, which is what we ul¬ 
timately care about? Obviously the difference will be 
negligible in stellar populations older than ^ 4 Myr, 
when the ionizing luminosity drops precipitously. Even 
in younger populations, the nebular contribution to the 
total luminosity is only ^ 10% for E275W and E438W fil¬ 
ters, and thus a difference in how this emission is treated 
is again unimportant. The nebular contribution is most 
significant for E336W, E555W, and E814W, where it can 
reach ^ 60 — 70% of the total (e.g., Reines et al. 2010). 

In the case of E555W, the nebular contribution is dom¬ 
inated by the [O ill] AA4959, 5007 and lines, with ad¬ 
ditional contributions from bound-free and two-photon 
continuum emission, and from the Ha and [N ll] A6584 
lines (which are intrinsically bright but lie near the edge 
of the filter response curve, and will be shifted out of 
the filter entirely for even modest redshifts). Over the 
observationally-plausible range of ionization parameters 
(roughly log 7/ = —3 to —2 - e.g., see the compilation 
in Verdolini et al. 2013), both Balmer line and bound- 
free production per ionizing photon injected vary by a 
factor of ^ 2. The [O ill] emission, as well as the other 
metal lines, are sensitive to both the temperature and the 
ionization state of the nebula. Consequently, they can 
easily vary by an order of magnitude as the ionization 
parameter changes (e.g., Kewley et al. 2001; Yeh et al. 
2013; Verdolini et al. 2013), and at the highest ioniza¬ 
tion parameters they can be a factor of several brighter 
than the HP line. Eor E814W, the nebular contribu¬ 
tion is dominated by free-free and bound-free continuum, 
with an a sub-dominant contribution of Paschen lines and 
[S III] AA9069, 9532. The two continuum sources and the 
Paschen lines are relatively insensitive to the ionization 
parameter, varying by less than a factor of 2. In contrast, 
the [S III] lines, while they are sub dominant, have the 
largest potential variation. Their strength can change by 
at least an order of magnitude over the plausible range of 
ionization parameters. Einally, the nebular contribution 
E336W is strongly dominated by the bound-free and two- 
photon emission. The former varies by less than a factor 
of 2 over the plausible ionization parameter range, while 
the latter stays almost entirely constant until the den¬ 
sity reaches ^ 10^ cm“^, at which point the two-photon 
process is collisionally quenched. 

Taken together, these factors suggest that differ¬ 
ences in the assumed ionization parameter between 
Yggdrasil and cluster_slug will lead to at most factor 
of ^ 3 differences in E336W, E555W, and E814W lumi¬ 
nosity in young stellar populations. These variations are 
not trivial, but we shall see below that they constitute 


a fairly minor uncertainty in comparison to those asso¬ 
ciated with degeneracies between age and extinction, or 
variations due to stochasticity. However, as a final cau¬ 
tion, we note that the prescriptions for nebular emission 
in both cluster_slug and Yggdrasil are only rough ap¬ 
proximations to the real complexity of H ll regions. Eor 
example, both cluster_slug and Yggdrasil assume a 
fixed, age-independent covering fraction. However, in 
reality H ll regions expand with time while the obser¬ 
vational aperture remains fixed, so the covering fraction 
should drop with age (cf. Whitmore et al. 2011), at a 
rate that depends on both the cluster’s luminosity and 
the density of its surrounding, both of which affect the 
rate at which its H ll regions expands. Thus age esti¬ 
mates below ^ 5 Myr should be regarded with extreme 
caution regardless of the fitting method. 

3. RESULTS FOR INDIVIDUAL CLUSTERS 

Throughout this section we will refer to the re¬ 
sults of an analysis of the cluster populations of 
NGC7793e, NGG7793w, and NGG628e performed us¬ 
ing cluster_slug. Eor the convenience of other authors 
who wish to use our results without having to rerun the 
full analysis, we summarize the cluster_slug output in 
Table 3. A full machine-readable version of this Table 
is included in the electronic edition of this paper, along 
with an analogous table containing the same results for 
the objects we have classified as unlikely to be genuine 
star clusters. We provide the latter as a service for those 
who might wish to make their own classifications. 

3.1. Are the cluster_slug Libraries Consistent with 
the Observations? 

As a first, most basic question, we examine the ex¬ 
tent to which the libraries of stochastic models repro¬ 
duce the colors and magnitudes found in the observed 
sample. The kernel density estimation method used by 
cluster_slug will return results even if the correspon¬ 
dence between the model libraries and the observations 
is very poor, but those results should be considered re¬ 
liable only to the extent that the models and observa¬ 
tions are in reasonable correspondence. To check whether 
this is the case, in Eigure 2 we show the distribution of 
our cluster_slug model libraries in color-color space, as 
compared to the observed catalog. As the plot shows, the 
model libraries generally cover loci in color-color space 
very similar to the observations. 

By itself, agreement in cuts in color-color space does 
not confirm that our libraries are a good analog to the 
observations. Even if the real and synthetic data show 
similar distributions in particular projections, they may 
be differently distributed in higher dimensions. More¬ 
over, unlike in deterministic models like Yggdrasil, the 
total mass and absolute magnitude are not free parame¬ 
ters in the cluster_slug models. Because the amount of 
stochasticity depends on how well the mass distribution 
is sampled, the dispersion in color at fixed age is a func¬ 
tion of the absolute magnitude. As a result, it is possible 
that our models cover the same range as the data in color- 
color space, but that they might not in color-magnitude 
space. 

To perform a more quantitative comparison of the ob¬ 
served and synthetic catalogs, we therefore examine the 
distribution of distances in photometric space between 
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Figure 2. Color-color plots comparing our model libraries 
against the observed star clusters in NGC7793e, NGC7793w, 
and NGG628. In all panels the y axis shows the color 
F336W — F275W, while the x axis shows the color indicated. 
(Note that, for NGG 628e we have used F435W in place of F438W 
- see Table 1.) Filled markers show every 20th cluster in the 
observed photometric catalogs, with the symbol shape and color 
indicating the galaxy as shown in the legend; the filters used for the 
photometric points are those indicated in Table 1. Golored points 
show 1% of the clusters from each of the cluster_slug libraries, 
selected at random; points shown for the cluster_slug libraries 
are all for WFG3 filters. Each row corresponds to one of the 
libraries listed in Table 2, as indicated in the right column. Points 
are colored by the age of the model, as indicated in the color bar. 


the observed star clusters and their closest analogs in 
our synthetic catalogs. We define the absolute and nor¬ 
malized photometric distances between an observed star 
cluster and the ith member of a cluster_slug library 

by 




2^norm_ 




-y' 

N ^ 


MFj,ohs — Mpj 


AM 


Fj ,obs 


(5) 

(6) 


where N is the number of filters, Mf-^ ohs is the magni¬ 
tude of the observed cluster in filter Fj, AMFj^obs is the 
observational error on this value, and Mf ,i is the mag¬ 
nitude of the ith library cluster in filter Fj. The exact 
filters used in this comparison are those given in Table 1, 
so generally N = 5; however, there are a very small num¬ 
ber of clusters which lie outside the image in one of the 
filters, and for these N = 4. 

For each observed cluster, we define and as 

the absolute and normalized distances to the 5th nearest 
neighbor in one of the dust er_slug libraries. (Note that 
the 5th nearest neighbor is not necessarily the same simu¬ 
lated cluster for the absolute and normalized distances.) 
Figure 3 shows the cumulative distribution of F 5 and 
2 ^norm each of our sample galaxies compared to each 
of our model catalogs; plots for the nth nearest neigh¬ 
bor, with n ^ 1 — 20, are qualitatively similar. Clearly 
nearly all of our observed star clusters have close analogs 
in our synthetic libraries. For our fiducial model library, 
pad_020_kroupa_MW, we find that 71% of observed clus¬ 
ters have at least 5 matches in our library within the la 
photometric errors (i.e., 71% have < 1), 95% have 

at least 5 matches within the 2a photometric errors, and 
99% have 5 matches within the 3 cr errors. The figures 
are comparable for the other model libraries, and the 
differences between them are quite minor. The normal¬ 
ized distances are generally smallest for NGC0628e and 
largest for NGG 7793e, but this is more a reflection of the 
size of the photometric errors in those fields than of any 
intrinsic differences between the cluster populations. 

While the nearness of our library points to the obser¬ 
vations is encouraging, we should inquire a bit further 
about what the distribution of nth nearest neighbor dis¬ 
tances should look like if our models are in fact a good 
fit. For a Poisson process, the cumulative distribution 
function of nth nearest neighbor distances r for a point 
at position x is given by 


n—1 


dn{u I CC) — 1 ^ ^ 

/c=0 


A(r, 

k\ 


(7) 


where A(r, x) is the expected number of points in a ball 
of radius r centered on x. Intuitively, this is simply the 
statement that the distribution of nth nearest neighbor 
distances is 1 minus the probability that there are n — 1 
or fewer points within a ball of size r around the point 
in question. Note that we are free to use any metric to 
measure r, so we can use the normalized photometric 
distance as well as the absolute one. Evaluating the ex¬ 
pectation value A requires knowledge of the shape of the 
underlying probability distribution around the point of 
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Figure 3. Cumulative distribution functions of photometric 
distances between the observed star clusters and the 5th nearest 
neighbor in our synthetic cluster_slug catalogs. The left column 
shows the raw 5th-nearest neighbor photometric distance 
(equation 5), while the right column shows i^norm (equation 6), 
the 5th nearest neighbor distance normalized to the photometric 
errors. Each row is for a different cluster.slug library, as 
indicated by the labels in the left column, and different line 
colors and styles are for different galaxies, as indicated by the 
legend. The black line, marked “All”, is the summed CDF for the 
three fields. The vertical dashed lines in the left column indicate 
D = 0.1 mag, the bandwidth we use for kernel density estimation. 

interest. In principle we could evaluate this for each of 
our points from the kernel density estimate for the PDF, 
but doing so would be quite computationally intensive, 
and for small distances would likely depend on our choice 
of bandwidth parameter. Instead, we can qualitatively 
check whether our nearest neighbor distribution is con¬ 
sistent with our models being a good match to the data 
by making a few simplifying assumptions that allow us 
to evaluate A analytically. 

Suppose that our library is large enough that we ex¬ 
pect there to be I simulations from the library within a 
normalized photometric distance = 1 of each ob¬ 

served point. Further suppose that the library of simula- 



Figure 4. Cumulative distribution function of 5th nearest 
neighbor photometric distances normalized by the observed 
photometric errors, The thick black line shows the 

distribution we measure summed over all our sample fields and 
using our fiducial model library, pad_020_kroupa_MW. The various 
color lines show the expected distribution of 5th nearest neighbor 
distances, computed using equations (7) and (8) for various values 
of the expected number of library models within the photometric 
error circle, and the dimensionality of the data, M, as indicated 
in the legend. 

tion points lies along an M-dimensional manifold within 
our 5-dimensional photometric space; for example, if the 
simulations near a particular point mostly lie along a 
line in photometric space, then the number within a 
given photometric distance increases linearly with dis¬ 
tance, and M = 1. If they are mostly along a plane then 
M = 2, and so forth. In this case we have 

X{r)=ir^, ( 8 ) 

and we can evaluate the expected nth nearest neighbor 
distribution d^(r) directly. Figure 4 shows how the mea¬ 
sured distribution of nearest neighbor distances for our 
fiducial model compares to the expected distributions 
with various plausible values of i and M. We see that 
the observed distribution is roughly consistent with the 
distribution we expect for ^ = 5 — 6, and M = 1, that is, 
the nearest neighbor distance distribution is about what 
we would expect if there were typically ^5 — 6 library 
models within the photometric error ellipse of each ob¬ 
servation, and if on small scales the distribution of points 
in photometric space mostly lies along a line. Thus our 
distribution of nearest neighbor distance is broadly con¬ 
sistent with what we would expect for a well-sampled 
library drawn from the same distribution as the data. 

To summarize, we find that the cluster_slug libraries 
are generally in excellent agreement with the observed 
photometric properties of the LEGUS sample fields. The 
majority of our sample has < 1, meaning that, for 

the typical observed cluster, there are at least 5 simulated 
clusters in the duster_slug libraries whose magnitudes 
in all filters are identical to within the photometric er¬ 
rors. We find with no obvious differences in the degree of 
agreement based on the choice of metallicity, evolution¬ 
ary tracks, IMF, extinction law, or nebular emission. 

3.2. Posterior PDFs of Star Cluster Mass, Age, and 
Extinetion 

Having verified that our synthetic libraries produce 
reasonable matches to the data, we now focus on our 
fiducial library, pad_020_kroupa_MW, and our fiducial 
prior distribution, pprior(logM, logT, Hy) oc , 












































Stochastic Star Cluster Models in LEGUS 


11 


leaving a discussion of the dependence of the results on 
these choices to the subsequent sections. We analyze 
the entire photometric catalog described in Section 2.1, 
and for each cluster we compute the marginal posterior 
probability of mass, age, and extinction on a grid of 128 
points each, covering the full range of each of these values 
present in our synthetic library. 

The computation is relatively fast - deriving each pos¬ 
terior PDF requires ^ 1 CPU-second per cluster for 
most high-confidence clusters in the catalog; those with 
the largest photometric error bars, or that are not well- 
fit by any models in our catalog (as is the case for 
many of the visually-rejected candidates, which we have 
nonetheless analyzed for completeness) may take up to 
a few tens of CPU-seconds, since large error bars re¬ 
quire that we search a larger volume of parameter space. 
Overall, we find that deriving posterior PDFs for all 
^ 3000 candidate clusters in the full catalog using a sin¬ 
gle cluster_slug library and sets of priors requires a 
few hours using a multi-core workstation; performing a 
similar analysis for the ^ 600 high-confidence clusters 
requires tens of minutes. 

Figures 5, 6, and 7 show sample marginal posterior 
distributions for mass, age, and extinction for 15 clus¬ 
ters per field in our 3 sample fields. The clusters shown 
are chosen to be uniformly distributed in the 50th per¬ 
centile estimates of their log mass, log age, and extinc¬ 
tion. From the plots, we can see that in most cases the 
cluster_slug models identify a fairly narrow range of 
possible masses for each cluster, with a typical interquar¬ 
tile range of ^ 0.2 — 0.3 dex on the posterior probabil¬ 
ity. The distributions are for the most part unimodal. 
However, we can see that there are a few cases where 
the posterior mass distribution is broader or even bi- 
modal, or where there is a tail of probability extending to 
very different masses, so that the 10th or 90th percentile 
whiskers extend very far beyond the 1st to 3rd quartile 
range. 

In comparison, the posterior probability distributions 
of age are somewhat broader and more likely to be bi- 
modal. In the bimodal cases there is often one peak at 
a relatively old age, and another at an age of ^ 10^’^ 
yr. This is particularly true for NGC628, which has the 
broadest photometric errors. The relative weighting of 
the two possible age fits, as we shall see below, is not in¬ 
dependent of our choice of priors. In these cases the me¬ 
dian age may not be a good representation of the actual 
age, because the median occurs near a local minimum of 
the PDF that lies partway between the two peaks. The 
posterior PDFs of extinction are also quite broad. In 
some cases they are bimodal, while in others they dis¬ 
play a single peak but with an extended tail. 

The nature of these bimodal fits is illustrated in Fig¬ 
ure 8, which shows various 2D projections of the posterior 
PDFs for one example cluster from NGC628e, which is 
typical of many of the bimodal fits we find. As the Fig¬ 
ure shows, the data are consistent with two “islands” 
of probability. The young island corresponds to an ex¬ 
tinction Ay ^ 0.2 — 1.5, age T < 3 — 10 Myr, and 
mass M ^ 3000 Mq, while the old one is centered near 
Ay - O.I, T - 500 Myr, M - 3-5 x 10^ M©. The color 
and luminosity of this cluster can therefore be fit well by 
either a relatively massive, extinction-free, old cluster, or 
a younger, somewhat less massive cluster that is red due 
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Figure 5. Box and whisker plot showing the marginal posterior 
probability distributions for star cluster mass for 15 example 
clusters per field in our 3 sample fields; the cluster IDs in the 
LEGUS photometric catalog are as indicated, and clusters are 
ordered from smallest to largest estimated 50th percentile mass. 
Note that ID numbers appear alternately above and below the 
boxes and whiskers. For each cluster, the blue colored band shows 
the probability density at each mass, as indicated in the color 
bar. Note that the color bar has been clipped above probability 
densities of 1.0 in order to reveal lower probability density features. 
The box and whisker plots (red) show percentiles: the lower and 
upper boxes indicate the range from the 1st to 2nd quartiles, and 
from the 2nd to 3rd quartiles, respectively. The lower and upper 
whiskers extend to the 10th and 90th percentiles, respectively. 


to greater extinction. 

The conclusions that many clusters when analyzed 
stochastically show multiple probability maxima is not 
new, and has been pointed out previously by Fouesneau 
et al. (2012, 2014), de Meulenaer et al. (2013, 2014, 2015), 
and Krumholz et al. (2015). Indeed, one could obtain 
such a result even with a deterministic method, provided 
that one used the full posterior PDF rather than an ap¬ 
proximation to it such as a Gaussian centered on the 
local minimum of The primary reason is that there 
are a number of places in color space where star clusters 
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Figure 6. Same as Figure 5, but now showing the marginal 
posterior distributions for age. Note that the clusters shown are 
not the same as the ones shown in Figure 5. 

with disparate physical properties are nonetheless very 
similar in color. The result is a likelihood function that 
is not well approximated by a uni-modal Gaussian. 

3.3. Dependence on Choice of Priors 

To what extent do the results for the posterior prob¬ 
abilities depend on the choice of prior probability dis¬ 
tribution? As noted above, the priors certainly matter 
for ages young enough that the color provides few con¬ 
straints, but the priors may also matter in other parts 
of parameter space as well. To answer this question, 
we continue to use our fiducial library, but now consider 
prior probability distributions of the form 

M^+^T, log(T/yr) < 6.5 . . 

M^+iT^+1, log(T/yr) > 6.5 ^ ^ 

with = — 1 or —2, and 7 = 0, —0.5, or —1; the combi¬ 
nation [3 = —2, 7 = —0.5 is our fiducial choice.That 

Note that the +l’s in the exponents in equation (9) occur 
because for cluster_slug we specify the priors on the log of mass 


Figure 7. Same as Figure 5, but now showing the marginal 
posterior distributions for visual extinction Ay- Note that the 
clusters shown are not the same as the ones shown in Figure 5. 

is, we consider cases where we take the prior on mass to 
be either fiat in log mass (/? = — 1 , all log masses equally 
likely) or with lower log masses more likely = — 2 , 
comparable to what is observed), and where we take the 
age distribution above ^ 3 Myr to be either flat in age 
(7 = 0 ), flat in log age (7 = — 1 ), or intermediate be¬ 
tween the two (7 = —0.5). Flat in age is what would be 
expected if clusters, once formed, never disrupt, while 
flat in log age is what would be expected if clusters dis¬ 
perse such that the survival probability is equal for each 
decade in time. We do not consider variations in the 
prior on Ay, as there seems to be little theoretical or 
observational motivation to do so. The effects of varying 
the prior on Ay are likely degenerate with the effects of 
varying the prior on age. 

In Figure 9 and Figure 10, we show how the masses 
and ages we derive for our star clusters depend on our 
choice of prior. We see that the choice of prior has rela- 

and age, while the conventional definitions of /3 and 7 are in terms 
of distributions of mass and age, rather than log of mass and age. 
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Figure 8. Triangle plot for an example cluster (ID 56) in 
NGC628e. Line plots show the marginal posterior PDFs for 
logM, logT, and Ay, while raster plus contour plots show 
the joint marginal posterior PDFs for the joint PDFs of these 
quantities in combination. Colors indicate probability densities as 
indicated in the color bar, and contours are spaced in intervals of 
0.2. All PDFs are normalized to have unit integral. 

lively little impact on the results for most cluster masses, 
typically moving the median by an amount significantly 
less than the 10th to 90th percentile range. The prior 
that deviates most from the others is/3 = — 2 , 7 = — 1 . 
With this choice the majority of clusters are still largely 
unchanged, but a substantial tail of clusters appears for 
which the deviation A logM is substantially less than 
zero, indicating that the prior /3 = — 2,7 = —1 leads 
us to a substantially smaller mass estimate than alter¬ 
native choices. The effect of prior choice on the inferred 
ages is somewhat greater, but as with masses, the only 
case that shows a very significant variation is P = — 2 , 
7 = — 1 , which again produces a tail of clusters with 
A log T substantially negative. In most cases this is sim¬ 
ply an understandable shift in power between two peaks 
of a bimodal PDF. As in the example shown above in 
Figure 8 , in some cases there are two islands of proba¬ 
bility that both fit the observed photometry reasonably 
well. In these cases a change in priors will enhance one 
island and suppress the other. Not surprisingly, in these 
cases the median shifts considerably. 

However, there are also a number of extreme outliers, 
where the amount by which the posterior PDF shifts 
when we change our priors is so great that the previ¬ 
ously inferred median now lies outside the lOth to 90th 
percentile range. Many of these points represent clusters 
with substantial photometric errors, clusters that are not 
well-fit by our model libraries (perhaps because they are 
composites of multiple populations of different ages), or 
both. It is not surprising that the derived results for 
these clusters should be very sensitive to the choice of 
prior. When the error bars on the observations are large, 
the likelihood ratio is close to flat, and thus the pos¬ 


terior probability distribution is little altered from the 
prior one. Something analogous occurs if no model clus¬ 
ter in our library is a good fit to the observations, and 
instead there are a broad range of models that are less 
accurate fits. 

However, there is also a population of clusters without 
large photometric error bars, and that are well-fit by our 
model library, that nevertheless show very substantial 
changes in their posterior PDFs depending on our choice 
of prior. To understand what is happening in these cases, 
in Figure II we show the 2D posterior PDF for an ex¬ 
ample cluster whose posterior PDF changes substantially 
depending on the choice of prior, but that does not have 
unusually large photometric errors, and that is well-fit by 
our models. As the plot shows, this cluster also has two 
islands of probability, one centered at a mass of ~ 10 ^-^ 
Mq and an age of ^ 500 Myr, and a second centered 
at a narrow range of ages ^5 — 10 Myr and masses of 
10^ — 10^ Mq. Unlike the degeneracy between age and 
extinction shown in Figure 8 , these two possibilities need 
not sit at substantially different Ay values. Instead, the 
degeneracy take a rather different form, which is more 
closely related to IMF sampling. 

The older age possibility corresponds to typical colors 
and luminosities for clusters of that mass and age range. 
On the other hand, the young age case corresponds to 
clusters that are on the far tail of the luminosity and color 
distributions for that age and mass. Thus the young age 
case represents extremely unusual photometric proper¬ 
ties for clusters so young and low mass to have, resulting 
from very improbable draws from the IMF. If all ages are 
considered equally likely, then the rare young, low-mass 
cases are rejected as unlikely. However, the difference 
between 7 = — I and 7 = 0, as indicated by red and blue 
colors in Figure II, amounts to the difference between a 
prior that says that there should be ^ 100 times as many 
clusters with ages from 5 — 10 Myr as clusters with ages 
of 505 — 510 Myr, and a prior that says there should be 
roughly equal numbers of clusters in those two ranges. If 
our prior is the former rather than the latter, then our 
Bayesian estimate assigns comparable total probabilities 
to the two possible fits, leading to a very different poste¬ 
rior PDF. 

It is an open question to what extent the dependence 
on the choice of prior could be reduced or removed by the 
availability of Ha data, which is particularly sensitive to 
the youngest ages and therefore good at discriminating 
between otherwise-degenerate models. Fouesneau et al. 
( 2012 ), analyzing star clusters in M83, found that Ha (as 
measured in their case by HST^s F657N filter) was very 
helpful in breaking degeneracies between fits. However, 
their analysis did not make use of F275W and instead 
had F336W as its bluest filter. The extra UV coverage 
provided by LEGUS’s F275W band should provide at 
least some of the same sensitivity to very young ages that 
Fouesneau et al. obtained from their Ha data. Moreover, 
Fouesneau et al.’s sample was limited to clusters with 
masses above ^ 10^ Mq, where stochastic effects are 
somewhat less important than in our sample. Since Ha 
is produced primarily by the most massive stars, it is 
particularly vulnerable to stochastic effects, which might 
reduce its ability to discriminate between models for our 
data set. In any case, a survey to obtain Ha data for a 
subset of LEGUS galaxies is underway (HST-GO-I3773, 
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Figure 9. Comparison of how the posterior masses we derive for the star clusters depend on our choice of priors. In each panel, the x axis 
shows the log of mass (logM)a; derived with one choice of prior (characterized by (/d, 7 ), the slopes of the mass and age distributions - see 
equation 9), while the y axis shows the difference A log M = (log M)y — (log M)x between this value and the 50th percentile derived using a 
different prior. In each column, all panels use the same prior to generate the value (logM)a; shown on the x axis, as indicated by the label 
at the top of the column; in each row, all panels use the same prior to generate (logM)^ and thus AlogM, as indicated by the labels at 
the right of the rows. The leftmost column, highlighted, uses our fiducial choice (^ = — 2,7 = 0) for (logM)^. Symbols of different shapes 
and colors correspond to different galaxies, as indicated in the legend, and the dashed horizontal line shows A log M = 0, indicating that 
the results are independent of the prior. Finally, the point plotted for each cluster is the inferred 50th percentile mass, while the horizontal 
and vertical error bars show the range from the lOth to 90th percentile as inferred for each set of priors. For the vertical error bars, the 
range in (AlogM) plotted is the range in (logM)y only, rather than reflecting a composite of (logM)^ and (logM)^. We show the lOth 
to 90th percentile range for only a subset of the data in order to minimize clutter. 
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Figure 10. Same as Figure 9, but showing a comparison of ages rather than masses derived using different priors. 
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Figure 11. Triangle plot for an example cluster (ID 320) in 
NGC628e. Panels are similar to those in Figure 8 . In the panels 
along the central diagonal showing the ID marginal PDFs, the 
black line corresponds to the fiducial case (/3 = — 2 , 7 = —0.5), 
while the blue and red lines refer to the alternate priors (3 = — 2 , 
7 = 0 and /! = —2, 7 = —I, respectively. In the three lower left 
panels, blue and red colors show the 2D PDF on a logarithmic 
intensity map (as indicated in the colorbar) for (3 = — 2 , 7 = 0 
and /5 = —2, 7 = —I, respectively. Black contours show the 2D 
PDF for the fiducial case, with contours placed at log probability 
density values starting at —I and increasing by 0.5 per contour. 
The ID and 2D PDFs are all normalized to have unit integral. 


PI: R. Chandar). As those data become available we will 
re-run our analysis pipeline using them, which should 
provide an answer to this question. 

It may also be possible to break degeneracies and re¬ 
move dependence on the choice of prior using other dis¬ 
criminators. For example, the amount of scatter in sur¬ 
face brightness within the aperture (Whitmore et al. 
2011 ) and the number of red stars in the vicinity of the 
cluster (Kim et al. 2012) have both been suggested as 
age-indicators. At present it is not clear how to build 
these into a Bayesian framework such as the one we have 
developed, as this would require extension of the likeli¬ 
hood function to include this information. 

3.4. Dependence on Choice of Tracks, IMF, Metallicity, 
Extinction Curve, and Nebular Emission 

We next examine the extent to which our results de¬ 
pend on our choice of tracks, IMF, metallicity, and 
extinction curve. In Figure 12 we show a compari¬ 
son between results derived using our fiducial model, 
pad_020_kroupa_MW, and results derived using models 
with different extinction curves, metallicities, IMFs, and 
stellar tracks, and omitting nebular emission. In all 
these comparisons we use our fiducial priors, /3 = — 2 , 
7 = —0.5, but the results are comparable for other priors 
provided that we use the same prior for each library. Ex¬ 
amining the figure, it is clear that the choice of IMF and 
extinction curve make almost no difference to the final 
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Figure 12. Comparison of posterior PDFs derived using our fidu¬ 
cial model, pad_020_kroupa_MW, and various alternatives: from 
top to bottom, pad_020_chabrier_MW (Chabrier (2005) IMF in¬ 
stead of Kroupa (2001) IMF), pad_020_kroupa_SB (starburst in¬ 
stead of Milky Way extinction curve), pad_020_kroupa_MW_noneb 
(same as the default library, but omitting nebular emission), 
gen_0I4_kroupa_MW (Geneva stellar tracks at Z = 0.014 versus 
Padova tracks at Z = 0 . 02 ), pad_004_kroupa_MW {Z = 0.004 in¬ 
stead oi Z = 0.02), and pad_008-kroupa_MW {Z = 0.008 instead 
oi Z = 0.02). In each panel, the fiducial model is on the x axis, and 
the comparison is on the y axis. The left column shows masses, 
and the right shows ages. Points plotted are the 50th percentile 
estimates of each quantity, and error bars (plotted for only a subset 
of the points to reduce clutter) indicate the 10th to 90th percentile 
range. Different colors and plot symbols correspond to different 
fields, as indicated in the legend. 
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posterior PDF. This is to be expected. For the IMF, the 
Chabrier (2005) and Kroupa (2001) IMFs we have tried 
differ mostly at the brown dwarf end, which has little im¬ 
pact on the integrated light properties even for old stellar 
populations. They also differ in that one is truncated at 
100 Mq and the other at 120 Mq, but those stars ap¬ 
pear to be rare enough that the difference made to the 
integrated light by including or omitting them is smaller 
than the level of variation induced by IMF sampling ef¬ 
fects. As for Ay, recall that the Ay values we infer are 
generally modest. If there is little extinction, then the 
shape of the extinction curve matters little. The choice 
of IMF also makes little difference, which is, again, not 
surprising. 

Whether we include nebular emission or not only 
makes a difference at young ages and low masses, though 
the latter is a selection effect - our ffux limit is such 
that our low-mass clusters are exclusively young. In this 
range, we find that models excluding nebular emission 
produce systematically higher masses. This result is easy 
to understand: if we ignore the light produced by the 
nebula, then a higher stellar mass is required to pro¬ 
duce the observed light. The effect of ignoring nebular 
emission on ages is more subtle, and again points to the 
importance of priors. We find that models excluding the 
nebular light do not produce 50th percentile ages below 
^ 3 Myr, while our fiducial models show no such ex¬ 
clusion. The result for the no-nebular case is driven by 
the fact that the colors of stars are essentially constant 
at ages < 3 Myr, so the likelihood function we compute 
by comparing to the observed colors is fiat over this age 
range. Combined our prior that all ages T below 10^'^ yr 
are equally likely, the posterior distribution with respect 
to logT must therefore peak at larger T. On the other 
hand, when including nebular emission colors do evolve, 
at least mildly, at younger ages, and thus the likelihood 
function is not fiat and we differentiate ages below 3 Myr. 
As a result, the 50th percentile values populate the full 
range of ages < 3 Myr. 

The choice of tracks has modest but non-negligible ef¬ 
fects. Compared to our fiducial case, the Geneva tracks 
seem to favor somewhat lower masses and ages. However, 
it is important to notice that the changes are largest for 
those clusters that have very significant uncertainties, so 
that the one-to-one line still falls mostly within the 10th 
to 90th percentile range. In effect, in those cases where 
the colors are ambiguous and the posterior PDF is multi- 
peaked, changing from one set of tracks to another can 
favor or disfavor one of the two peaks. 

The largest effect is from metallicity. Compared to 
Solar metallicity, the lower metallicity models favor sub¬ 
stantially lower ages and masses for intermediate mass 
and age clusters, and slightly higher masses and ages 
for the youngest and lowest mass clusters. The effect at 
young ages and low masses appears to be similar to that 
seen in the the no-nebular versus nebular comparison. 
This results from a complex series of effects of metallicity 
on the color: lower metallicity increases the ionizing ffux 
and raises the nebular temperature, thereby increasing 
continuum and hydrogen recombination line emission. 
On average it also weakens emission from metal lines, 
though some lines that are particularly temperature- 
sensitive (e.g., [O III] A5007) may also get brighter. The 
resulting shift in colors appears to favor older ages and 


higher masses for the youngest and lowest mass clusters. 
The effect at higher mass and age comes largely from the 
interaction of colors with IMF sampling effects. Choos¬ 
ing a low metallicity tends to drive the model colors to 
the blue, beyond the point where they agree well with 
the observations. To match the observed colors requires 
driving the colors back to the red, and one way of doing 
this is to have a lower mass cluster that under-samples 
the massive end of the IMF, thereby producing redder 
colors. 

Beyond the individual model comparisons, perhaps the 
most striking result shown in Figure 12 is how little dif¬ 
ference the choice of models makes. The only one of the 
variations that we have examined that appears to make 
a substantial difference is using Z = 0.004, and, at the 
youngest ages perhaps, including or not including nebu¬ 
lar emission. The primary reason for this insensitivity is 
that widths of the posterior PDFs for age and mass esti¬ 
mates are sufficiently large that they swamp systematics 
such as the choice of extinction curve, IMF, evolution¬ 
ary tracks, and, at least over a limited range, metallicity. 
Perhaps these choices would become important if we had 
additional observational constraints beyond the five pho¬ 
tometric bands to which we currently have access, but at 
present they are not the limiting factors in our knowl¬ 
edge. 

3.5. Comparison of Stochastic and Non-Stochastic 
Models 

Our final analysis step is a comparison of the re¬ 
sults we obtain using cluster_slug stochastic models 
to those we get from the deterministic Yggdrasil code. 
For this comparison, we use our fiducial library, 
pad_020_kroupa_MW, and prior distribution, (3 = —2, 
7 = —0.5. We can investigate this question on two 
levels. We can first ask how well cluster_slug and 
Yggdrasil models match the observed photometry, and 
whether one provides better matches of the observations 
than the other. We can then ask how the results we 
derive for cluster properties, and the nominal errors on 
those results, compare for the two methods. 

3.5.1. How Well Do cluster_slug and 
Yggdrasil Reproduce the Observed Photometry? 

To investigate the consistency or lack thereof between 
the cluster_slug and Yggdrasil models and the ob¬ 
served photometry, we must first develop a statistic to 
quantify it. For duster_slug, we parameterize the 
quality of agreement between the observations and the li¬ 
brary using normalized photometric distance 

between an observed cluster and its 5th closest neighbor 
in the cluster_slug library (equation 6), exactly as in 
Section 3.1. For Yggdrasil, the best fitting model is 
determined via a minimization procedure. From the 
value of determined for the best fit (and the number 
of degrees of freedom in the model), we can compute the 
Q statistic Q(x^), which formally is the probability that 
the value of would exceed the measured one even if 
the model were correct. Thus for example Q(x^) = 0.3 
means that, even if the best fit model were true, there is a 
30% chance that each observation of that cluster would 
yield photometry for which the x^ value would exceed 
our measured value. Thus values of Q(x^) near unity 


18 


Krumholz et al. 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0.0 0.1 0.2 0.3 


iir- (cluster_slug) 


Figure 13. Distribution of errors between the observed photom¬ 
etry and the models. In the central panel, for each cluster the 
X coordinate shows Dgorm (equation 6), while the y coordinate 
shows Q(x^). Values of Q < 10“^ are shown as upper limits, 
indicated by arrows. The histograms in the flanking panels show 
the binned distribution of Dnorm (^^p p^nel) and Q (right panel) 
values for all clusters; the lowest, hatched bar in the frequency 
distribution for Q shows the frequency of points with Q < 10 “^. 
In both histograms, the solid lines indicate the distribution of 
values we would expect for cluster.slug or Yggdrasil models 
that provide good fits to the data - see main text for details. 


correspond to good fits, while ones with Q(x^) ^ 1 indi¬ 
cate either that the error bars have been underestimated 
or that the best fitting model in the Yggdrasil library 
is probably not a good fit to the data. 

Figure 13 shows the distribution of Q(x^) versus 
for our photometric catalog. In this plot, good fits 
for cluster_slug lie on the left side of the diagram, 
while poor fits lie on the right. Similarly, good fits 
for Yggdrasil lie at the top of the plot, and poor fits 
toward the bottom. The plot shows, consistent with 
Figure 3, that data are generally well-reproduced by 
the cluster_slug library. The solid line in the his¬ 
togram of H^orm j^odels shows the expected distribu¬ 
tion for a library with ^ = 6 models within an obser¬ 
vational error circle and M = 1 dimensional data (see 
Section 3.1); it is simply the binned, differential version 
of the cumulative distribution shown in Figure 3. The 
cluster_slug distribution of errors is roughly consistent 
with this distribution, confirming that most clusters have 
good cluster_slug fits. 

The same is not true of the Yggdrasil models. For a 
good fit the Q(x^) statistic should be distributed so that 
a fraction / of models have Q(x^) < /• fact. Figure 13 
shows that there is a significant excess models at Q ^ 1. 
The solid line in the histogram of Q values shows the 
distribution we would expect for a good fit, and there 
are clearly more small values of Q than expected, and 
significantly fewer values close to unity. In particular, 
nearly 10% of models have Q < 10“^, while for a good 
fit such large discrepancy should occur 0.1% of the time. 
In contrast, < 5% of clusters have > 2, and only 

1% have > 3. 

One possible explanation for the excess of clusters with 
Q ^ 1 is that our observed catalog includes some objects 


that are not truly simple stellar populations, for example 
because they are blends of two clusters of different ages. 
However, we can immediately reject this as the domi¬ 
nant explanation, because such clusters should also be 
fit poorly by cluster_slug models, placing them in the 
lower right corner of the main panel of Figure 13. While 
there are a few objects of this sort, there are also a very 
large number of clusters for which Yggdrasil returns a 
Q value of ^ 10“^ or less, but cluster_slug neverthe¬ 
less finds at least 5 simulated clusters whose photom¬ 
etry matches the observed values within 1 or 2 times 
the photometric errors. Indeed, even for the subset of 
clusters for which Yggdrasil produces Q < 10“^, the 
mean value of unorm This is despite the fact that 

Yggdrasil is only fitting the colors (since the mass and 
thus the absolute luminosity are left as free parameters), 
while cluster_slug is fitting both the colors and the ab¬ 
solute magnitude. In contrast, there is no corresponding 
population of clusters in the upper right part of the dia¬ 
gram, which would indicate that Yggdrasil finds a good 
fit but cluster_slug does not. 

Many of the clusters for which Yggdrasil does not 
find good fits are those for which cluster_slug assigns 
relatively low masses; the mean of cluster_slug’s 50th 
percentile mass estimate for the subset of clusters for 
which Yggdrasil finds Q < 10“^ is 1000 Mq, smaller 
than the 50th percentile mass for the entire sample by 
a factor of several. This strongly suggests that much 
of the problem is driven by Yggdrasil’s assumption of 
a well-sampled IMF, which fails in the low mass clus¬ 
ters. For these cases, exploring the full range of stochas¬ 
tic color variation is crucial to finding a good match to 
the data. However, there also some relatively massive 
clusters for which cluster_slug finds a good match but 
Yggdrasil does not. These may be cases where, despite 
the cluster’s overall relatively large mass, its colors are 
still significantly influenced by a small number of stars 
undergoing short-lived phases of stellar evolution. For 
example, a small number of He core stars undergoing 
blue loops might dominate the UV light budget of an 
otherwise old and red stellar population. For these rare 
phases of stellar evolution, the continuous assumption 
used in Yggdrasil may fail even for relatively massive 
clusters. 

3 . 5 . 2 . Comparison of cluster_slug and Yggdrasil Best 
Fits and Errors 

Having discussed how well the models libraries match 
the data, we now ask how well the predicted cluster 
properties and the errors on those properties agree be¬ 
tween the two methods. Figure 14 shows a compari¬ 
son of the best-fit Yggdrasil values to the 50th per¬ 
centile cluster_slug values. This is not quite a fair 
comparison, since in some cases the 50th percentile 
cluster_slug mass does not actually lie near a prob¬ 
ability maximum. Nonetheless, the obvious alternative, 
plotting the cluster_slug point at the probability max¬ 
imum, is no better. Some imprecision is inevitable when 
comparing a full posterior PDF to a single best fit. 

Comparing masses in Figure 14, we notice that the 
models generally agree fairly well at high masses, but 
at low masses the Yggdrasil models on average tend 
to produce lower masses compared to cluster_slug’s 
50th percentile, making the scatter about the 1 — 1 line 
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Figure 14. Comparison of cluster masses (top panel), ages 
(middle panel), and extinctions (bottom panel) computed with 
Yggdrasil and cluster_slug. In all panels, the central point is 
plotted at the 50th percentile value returned by cluster_slug on 
the X axis, and the best fit returned by Yggdrasil on the y axis. 
Error bars, which we show only on a subset of points for clar¬ 
ity, indicate the 68% confidence interval for Yggdrasil, and the 
16th to 84th percentile range for cluster_slug. Point colors and 
shapes indicate the field for each cluster, and black dashed lines 
show the 1 — 1 relation. In the top two panels, large yellow points 
show the mean logarithmic mass and age, respectively, for both 
cluster.slug and Yggdrasil. Mean masses are computed by av¬ 
eraging the clusters in bins of 50th percentile cluster_slug mass 
from log(M/MQ) = 2 — 3, 3 — 4, 4 — 5, 5 — 6, and 6 — 7; ages are com¬ 
puted via an analogous procedure over bins from log(T/yr) = 6 — 7, 
7 — 8, 8 — 9, and 9 — 10. In the middle panel, the alignments of the 
Yggdrasil fits at a set of common ages is an artifact of the fitting 
procedure. 


somewhat asymmetry. Despite this, the mean masses 
for the population as a whole (indicated by the yel¬ 
low circles in Figure 14) are quite similar. The origin 
of this effect is almost certainly IMF under-sampling. 
For ^ 500 Mq clusters, the most common outcome by 
number is that the cluster will under-sample the mas¬ 
sive end of the IMF, and thus will be under-luminous 
compared to what would be expected for a fully-sampled 
IMF. Because Yggdrasil assumes full sampling it as¬ 
sumes a mass to light ratio that is too low, and ends up 
with a mass estimate that is too low as well. In contrast, 
cluster_slug correctly accounts for the sampling effect 
and assigns a broad PDF whose median value lies at 
higher masses than Yggdrasil’s best fit. The effect be¬ 
gins to appear at masses below roughly 10 ^-^ Mq, consis¬ 
tent with earlier analyses (e.g., Elmegreen 2000; Cerviho 
& Luridiana 2004, 2006; da Silva et al. 2012; Fouesneau 
et al. 2014). However, when averaging over a large pop¬ 
ulation of clusters, there will be a few that happen to 
be over-luminous rather than under-luminous for their 
mass, so that the mean is the same as for fully-sampled 
IMF. Consequently, while there is an asymmetric bias in 
the ages of individual clusters, estimates for the mean of 
the entire population as much less affected. 

There are also systematic offsets in Yggdrasil and 
cluster_slug’s age and Ay distributions. Here the phe¬ 
nomenology is more complex. As we have seen, the 
ages one derives from cluster_slug are not indepen¬ 
dent of the choice of priors. There are often genuine 
ambiguities, with multiple age-Ay combinations provid¬ 
ing plausible fits to the data. The 50th percentile value 
that cluster_slug derives will depend on how the priors 
weight these two reasonably good fits. The best-fit pro¬ 
cedure used by Yggdrasil should be roughly equivalent 
to searching for the highest peak in age-Ay space, based 
on priors that are fiat in logT and Ay. The default 
cluster_slug priors used for the results shown in Fig¬ 
ure 14 have priors that are flat in age rather than log age 
(though this is partly offset by the prior to low masses, 
which implicitly favors younger ages), and so it is not 
surprising that they produce somewhat different results. 
Despite the offsets, however, we see that the population 
means are again fairly similar between cluster_slug and 
Yggdrasil. 

The greatest discrepancy would appear to be in the 
values of Ay derived by the two methods. However, this 
is somewhat misleading, as the error bars on Ay are often 
extremely large. Thus despite the seeming large offset 
between the cluster_slug medians and Yggdrasil best 
fits, in reality almost all of the data points have the I — I 
lines within their allowed error range. 

Finally, we note that the uncertainty range produced 
by cluster_slug is in almost all cases larger than the 
one produced by Yggdrasil. This is particularly true 
of masses at the low mass end, where cluster_slug’s 
68% confidence range is often close to an order of 
magnitude wide, while Yggdrasil’s is so small that 
it is nearly invisible in Figure 14. This is likely be¬ 
cause Yggdrasil’s error estimate only includes the for¬ 
mal error coming from propagation of the photomet¬ 
ric errors though the deterministic model grid, while 
cluster_slug properly captures the additional uncer¬ 
tainties coming from stochasticity and from degenera¬ 
cies. Clearly the latter two effects dominate the total 
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error budget, leading cluster_slug to have much larger 
output errors than Yggdrasil. 

4. PROPERTIES OF THE CLUSTER POPULATIONS 

Having analyzed the properties of the individual clus¬ 
ters and their statistical and systematic uncertainties, we 
now seek to extent our analysis to the properties of the 
cluster population as a whole. Before commencing this 
exercise, we note that the results we derive in this section 
should not be taken as representative of the underlying 
properties of star clusters in our sample galaxies. In¬ 
stead, we seek to derive the properties of those clusters 
that have made it into the LEGUS photometric cata¬ 
log, which is the convolution of the true distribution of 
properties with the catalog selection. To convert these 
to intrinsic properties our analysis would have to be cor¬ 
rected for completeness, and such a correction is beyond 
the scope of the present work. In general we note that, 
because our catalog consists of visually-inspected clus¬ 
ters, and the main criterion for being subject to visual 
inspection is exceeding a magnitude limit, our sample 
should be closer to magnitude-limited than mass-limited. 
This will tend to result in a steeper age distribution than 
is present in reality (e.g.. Earners 2009). 

4.1. Methods 

Our goal here is to derive the probability distribution 
functions for ages and masses of clusters, or their joint 
distribution, given the posterior PDFs for each individual 
cluster that we have derived in the previous section. De¬ 
riving this requires some care. The traditional method 
to derive the properties of the population, when each 
cluster has been fit using a traditional approach, is to 
create bins in mass and age, and then assign clusters to 
the bin corresponding to the best fit. The errors on each 
bin can then be computed from Poisson statistics.This 
is the method we will use to analyze the Yggdrasil re¬ 
sults below. As usual, one must make a choice of how 
to place the bins, and we will consider two cases: one 
where the bins are placed uniformly in the log of mass or 
age, and one where they are placed so as to have a fixed 
number of clusters per bin. Note that the latter choice 
also maintains a fixed value of the relative error on the 
PDF: since the error for a bin containing N clusters is 
a/]V, the ratio of the error to mean is just 1/V^, and 
thus a fixed value of N implies a fixed fractional error. 

This method is clearly not well-suited to our posterior 
PDFs, which are broad and multiply-peaked. Moreover, 
we wish to avoid binning, since binning necessarily en¬ 
tails the loss of information. Instead, we proceed via a 
bootstrap resampling procedure. Consider the case of 
the mass distribution; distributions of age, extinction, or 
any combination of the three can be treated in an anal¬ 
ogous fashion. Our central estimate for the mass distri¬ 
bution of the population as a whole is simply the sum of 
the posterior PDFs of mass for all the observed clusters. 
That is, in a galaxy containing N clusters, for which we 
have determined a posterior probability distribution for 

A slightly more accurate procedure is to distribute the clusters 
between bins based on the convolution of the bin with their formal 
uncertainties, but since we have already seen that traditional 
method greatly underestimates the uncertainties (Section 3.5), this 
is only a marginal improvement. 


mass Pi{m) for the ith cluster, our central estimate of 
the population PDF is simply 

1 ^ 

Ppop(m) = — J2pi{m). (10) 

To estimate the uncertainty on this estimate, we per¬ 
form Nt trials. For each trial we draw N clusters from 
our observed sample with replacement, meaning that we 
may draw the same cluster more than once. We then 
compute the population PDF Ppop,t(^) for that trial via 
equation 10. Thus after Nt trials we have Nt population 
PDFs, and at each mass m we have Nt estimates for the 
value of Ppop,t(m). We can then compute a 90% confi¬ 
dence interval on Ppop(m) as simply the range from the 
5th to the 95th percentile of these Nt estimates, and sim¬ 
ilarly for any other confidence interval of interest to us. 
The result of this procedure is therefore both a central 
estimate and a confidence interval on Ppop(m) at each 
mass m, obtained without the need for binning, and us¬ 
ing the full posterior PDFs from our Bayesian analysis 
of the individual clusters. 

4.2. Marginal Mass and Age Distributions 

Figure 15 shows the inferred marginal mass distribu¬ 
tions for the entire population of star clusters, summing 
over all three of our fields. We emphasize that these plots 
show p(log M) rather than p{M), so that a flat line corre¬ 
sponds to equal numbers of cluster per logarithmic mass 
interval. We see that, for our cluster_slug based mod¬ 
eling, the mass distribution is consistent with a powerlaw 
dNjdM or perhaps very slightly steeper, over 

the mass range from roughly lO^’^ —10^’^ Mq. This result 
is essentially independent of the choice of model library 
or prior, and with relatively small statistical uncertainty 
based on bootstrap resampling. At lower masses the 
mass distribution flattens, almost certainly as a result of 
incompleteness, and at the lowest masses, M < 10^ Mq, 
the inferred posterior distribution depends fairly strongly 
on the prior. This is a result of posterior PDFs for indi¬ 
vidual clusters in this mass range being quite broad, so 
that the total shape depends on the prior. In compar¬ 
ison, the choice of library matters little at low masses. 
Our central estimate for the shape of the mass function 
at the very high mass end, approaching 10 ^ Mq, depends 
strongly on both the choice of library and the priors, but 
the statistical uncertainty is so large (due to the small 
number of clusters) that the bootstrap confidence inter¬ 
vals are mostly overlapping. 

Figure 16 shows the corresponding marginal age dis¬ 
tributions. As with Figure 15, we are plotting p(logT) 
rather thanp(T) (in contrast, for example, to the similar 
Figure 14 of Fouesneau et al. 2014), so that a flat line 
indicates a distribution where the cluster survival prob¬ 
ability per decade in age is constant. Absence of cluster 
disruption would appear as a line of slope +1 on this 
plot. Based on the cluster_slug analysis, these are dis¬ 
tributed roughly as dN/dT ^ T~^ at ages from roughly 
10^ — 10^’^ yr, with some dependence of the slope on the 
choice of prior and on the model library. The largest 
outlier is the result derived from the Geneva library, but 
this should probably be regarded as less reliable in the 
relevant age range due to its omission of TP-AGB stars. 
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Figure 15. Inferred posterior PDF for the mass distribution of the 
entire population of star clusters, marginalized over age and extinc¬ 
tion, and summed over the three sample fields. In the top panel, 
the different colored lines indicate central estimates for this quan¬ 
tity derived from cluster.slug using the different model libraries: 
the fiducial one (pad_020_kroupa_MW), one using a Chabrier 
instead of Kroupa IMF (pad_020_chabrier_MW), one using a star- 
burst instead of Milky Way extinction curve (pad_ 020 _kroupa_SB), 
one omitting nebular emission (pad_020_kroupa_MW_noneb), one 
using Geneva instead of Padova tracks (gen_020_kroupa_MW), 
and two with metallicities of Z = 0.004 and Z = 0.008 instead of 
Z = 0.020 (pad_004_kroupa_MW and pad_008_kroupa_MW). In 
the middle panel, we show central estimates for the mass PDF 
estimated using our fiducial library (pad_020_kroupa_MW), but 
a number of different priors (/ 3 , 7 ), as indicated. In the bottom 
panel, we repeat some of the models from the top two panels, this 
time with shaded regions indicating the 5th to 95th percentile 
confidence interval based on bootstrap resampling with 50,000 
trials. Also in the bottom panel, circles with error bars show 
the values inferred by binning the best fit masses inferred by 
Yggdrasil, using either 25 adaptive bins with equal numbers of 
clusters per bin (white points), or uniform bins spaced at intervals 
of 0.5 dex (salmon points). Error bars on the Yggdrasil points 
show the Poisson error. The black dashed and dotted lines that 
appear in all panels show the priors in mass corresponding to 
/3 = —2 and /3 = — I. 



log T [yr] 


Figure 16. Same as Figure 15, but now showing the age 
distribution instead of the mass distribution. Black lines show 
different priors in age corresponding to 7 = 0, —0.5, and —1.0 at 
ages above 10 ®'^ yr. 


The results with all models and libraries also show an 
excess of clusters at ages around 10^’^ yr. This is partly 
driven by the choice of prior, which is flat in age T (and 
thus rising in logT) up to this age. We have argued that 
this is in fact the correct prior to use in this age range, 
since cluster disruption will be negligible for clusters this 
young. However, the bump at 10^’^ yr might also be af¬ 
fected by completeness, which will preferentially suppress 
clusters at ages older than this due to the rapid decline 
in luminosity at older ages. At ages much below 10^’^ 
yr, the results appear to be quite sensitive to the choice 
of metallicity and tracks, and the treatment of nebular 
emission. The feature driving this is almost certainly the 
variations in nebular emission, which is affected directly 
by metallicity, and indirectly (though the ionizing lumi¬ 
nosity) by the choice of tracks. Eventually all models 
converge to our prior, dN/dT ^ constant, because, at 
ages below 1 Myr, even if we include nebular emission 
then there is no change in cluster colors with age. Since 
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the colors do not distinguish between ages < 1 Myr, the 
posterior at such young ages simply reduces to the prior. 

We also see a drop off in the cluster population at ages 
above ^ 10 ^-^ yr. At least some of this drop must oc¬ 
cur because these clusters’ low luminosities place them 
below our sensitivity limit, but there may also be a real 
decline in cluster numbers at older ages as well. If so, this 
would imply an increase in the rate of cluster disruption 
at ages around 10^-^ yr. However, given our incomplete¬ 
ness, it would be premature to draw any conclusions at 
this point. The shape of this drop off depends some¬ 
what sensitively on the choice of prior, with different but 
still physically-plausible priors giving shapes in this age 
range that can be outside each others’ 90% confidence 
intervals. It does not depend strongly on the choice of 
library, again with the exception of the Geneva models, 
which produce quite a different result but are likely less 
reliable. On the other hand, it is possible that the Padova 
libraries may also have difficulties in this age range, since 
it is particularly susceptible to choices in how one treats 
the hard-to-model TP-AGB phase. 

As noted above, to turn the distributions shown in Fig¬ 
ure 15 and Figure 16 into mass and age distributions for 
clusters in general will require completeness correction, 
a topic that we leave to a future paper. Moreover, the 
completeness correction will quite different for NGG628 
than for NGC 7793 due to their different distances 10 
Mpc versus ^ 3.5 Mpc). Our goal is therefore not to 
derive the properties of the cluster populations in these 
two galaxies, but rather to understand how such deter¬ 
minations are infiuenced by the method used to convert 
the photometry into physical properties. We can already 
see from the results obtained thus far that the age dis¬ 
tribution, at least in certain age ranges, will not be in¬ 
dependent of the choice of prior, or the choice of library. 

Comparing the cluster_slug results to the ones ob¬ 
tained using Yggdrasil, we find that the mass PDFs 
are roughly consistent within the error bars, though this 
depends to some extent on how one chooses to bin the 
Yggdrasil results, and on the prior that one chooses for 
the cluster_slug analysis. At small masses where the 
statistical uncertainties are largest, the Yggdrasil re¬ 
sults appear to be closest to what one obtains using 
cluster_slug with a prior distribution (3 = — 2 , 7 = — 1 . 
For the age distributions there is more discrepancy be¬ 
tween the Yggdrasil and cluster_slug results. The 
Yggdrasil data show sharp drops in the number of clus¬ 
ters at around and yr, with excesses near 

10 ^’^ and yr, though the prominence of these fea¬ 
tures depends on how one chooses to bin. The same 
structure in the Yggdrasil-determined ages is visible as 
the horizontal banding in Figure 14, and appears to re¬ 
sult from the choice to assign a single best fit even in 
parts of color space where colors provide relatively poor 
constraints on the true age. The cluster_slug mod¬ 
els correctly capture the uncertainty in such regions by 
spreading out the posterior probability, and show no 
corresponding features. Instead, the population PDF 
is smooth. Nonetheless, the broad distribution of ages 
in the Yggdrasil models is not dissimilar from that 
obtained with cluster_slug. Apparently while the 
stochastic and deterministic methods produce relatively 
large differences on a cluster-by-cluster basis, they agree 
much more closely when analyzing the cluster population 



Figure 17. Joint posterior distributions of age and mass for all 
clusters in all fields, marginalized over extinction. Colors and 
solid contours show the probability density in the (log M, log T) 
plane; PDFs are properly normalized, so the integral of the PDF 
over the plane is unity. The top six panels show the results 
for cluster_slug using different priors (/3, 7 ) as indicated, while 
the bottom panel shows the results for Yggdrasil, which have 
been determined by taking the Yggdrasil best fits and making 
a 2D histogram using bins 0.25 dex wide in both the mass 
and age directions. Dashed contours show the results using 
cluster_slug with our fiducial prior (/3 = —2, 7 = —0.5), and are 
the same in every panel in order to facilitate comparison. Both 
solid and dashed contours start at a log probability density of —3, 
and are spaced at intervals of 0.5 thereafter. 

as a whole. 

4.3. Joint Mass-Age Distribution 

In addition to examining the mass and age distribu¬ 
tions separately, we can examine them jointly. Figure 17 
shows the joint distributions of mass and age summed 
over all three fields, marginalizing over the extinction. 
All the panels shown use our fiducial model, but we ex¬ 
amine a range of priors. We also show the Yggdrasil re¬ 
sults, binned in the (log M, log T) plane. 

Qualitatively, the cluster_slug results with all priors 
and the Yggdrasil results agree that the observed clus- 
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ter population mostly lies along a band that extends from 
masses of ^ 10 ^ Mq and ages of ^ 10 ^-^ yr up to masses 
of ^ 10^-^ Mq and ages of ^ 10^-^ yr. This distribution 
is doubtless heavily influenced by the completeness of the 
survey, which would prevent us from detecting clusters 
that lie in the upper left of the diagram. On the other 
hand, the absence of clusters in the lower right corner, 
corresponding to massive, young systems, is a real fea¬ 
ture of the distribution, as such clusters should be readily 
detectable. 

Examining the results a bit more closely, we can notice 
that the choice of priors does influence the final distribu¬ 
tion in non-negligible ways. As the priors on mass and 
age become flatter, ^ = — 1 versus $ = —2^ and 7 = 0 
versus 7 = — 1 , clusters tend to shift from the lower left 
end of the populated band, at low masses and young 
ages, upward to higher masses and ages. Thus priors 
that favor low mass and young age, ;d = — 2,7 = — 1 , 
produce a prominent bump at masses of Mq 

and ages near 10 ^-^ yr, while flatter priors, P = — 1 , 
7 = 0 , remove this feature and distribute the clusters 
more evenly across a range of masses and ages. Elatter 
priors also produce a small island of probability at a mass 
of ^ 10^’^ Mq and an age of ^ 10^’^ yr. This island is 
produced by a handful clusters that are not well fit by 
either the cluster_slug or Yggdrasil models, making 
the cluster_slug results for them quite sensitive to the 
prior. 

The other very noticeable difference is between the 
Yggdrasil and duster_slug models. Here the key fea¬ 
ture is that the Yggdrasil-derived masses and ages are 
confined to a much smaller portion of the age-mass plane 
than the cluster_slug ones. This reflects the much 
smaller uncertainties that Yggdrasil assigns based on its 
fitting approach, as opposed to the Bayesian method 
of cluster_slug. This approach tends to concentrate 
the probability toward the peaks, suppressing the lower 
probability wings that are retained by using the full pos¬ 
terior. 

5. SUMMARY AND CONCLUSIONS 

We describe and investigate the performance of a 
Bayesian method to derive the masses, ages, and extinc¬ 
tions of star clusters observed using broadband photom¬ 
etry. Our sample data set for this analysis consists of 621 
visually-confirmed star clusters drawn from NGC628e 
and NGC 7793. These galaxies are part of the Legacy Ex- 
tragalactic UV Survey (LEGUS) sample, which provides 
HST photometry in the NUV, U, B, V, and I bands. Our 
method, implemented based on the slug software suite^^ 
and its cluster analysis tool dust er_slug, uses kernel 
density estimation coupled to implied conditional regres¬ 
sion to return the full marginal posterior probability dis¬ 
tributions of mass, age, and extinction. This technique 

For example, the cluster for which Yggdrasil assigns the 
highest mass, which appears in the right-most pixel in the 
Yggdrasil panel of Figure 17, has a Q parameter of 1.9 x 10“"^, 
indicating a very poor fit. The normalized photometric nearest- 
neighbor distance = 3.0 for our fiducial library, indicating 

that no cluster.slug models land closer than 3a to the cluster’s 
photometric properties; alternate libraries do no better. Given the 
strong disagreement with all sets of models, it seems likely that 
this object is not well described as a simple stellar population, or 
that all the model tracks we have available are deficient. 

https://www.slugsps.com 


proves to be particularly useful in two regimes. One is for 
low-mass clusters, below ^ 10 ^’^ Mq, where the initial 
mass function is not fully sampled and the relationship 
between physical and photometric properties is therefore 
non-deterministic, leading a broad posterior PDE that is 
not well characterized by a single best fit. The second 
case is for clusters that he at locations in color space 
where there are significant degeneracies between mass, 
age, and extinction, so that there are multiple possible 
fits of comparable likelihood at distinct sets of physical 
parameters. In both of these regimes, a full posterior 
PDE provides a more reliable result, and a more real¬ 
istic estimate of the errors, than methods that return 
a single best fit with a local error distribution around 
it. Indeed, by comparing our results to those produced 
by the non-stochastic cluster fitting code Yggdrasil, we 
find that the stochastic method removes some of the ar¬ 
tifacts found in the deterministic one, and that it returns 
substantially larger error distributions. 

At the level of individual clusters, our analysis tech¬ 
nique proves to be insensitive, within the plausible range 
of variations, to our choice of metallicity, extinction 
curve, stellar initial mass function, and evolutionary 
tracks. On the other hand, because in many cases the 
observed photometry is consistent with more than one 
combination of mass, age, and extinction, our choice of 
prior probabilities does affect the final result. This sen¬ 
sitivity occurs because, when more than one mass-age- 
extinction combination is consistent with the observed 
photometry, the relative weight we assign to different “is¬ 
lands” of probability depends on our priors. This means 
that, for an individual cluster, before drawing any firm 
conclusions about its age, or to a less extent its mass, one 
should be careful to ensure that the results are robust 
against variations in priors, whether explicit or implicit. 
It is possible that the addition of extra photometric data, 
particularly Ha, might break some of these degeneracies, 
but that remains to be seen. 

At the level of a cluster population, we find that 
the mass distribution in the range ^ 10 ^ — 10 ^ Mq is 
very robust against changes in either prior or assumed 
stellar models, and to whether we derive the masses 
using stochastic cluster_slug models or deterministic 
Yggdrasil ones. At lower masses, priors begin to mat¬ 
ter because the large amount of stochastic variation in¬ 
duced by incomplete IME sampling means that photome¬ 
try no longer provides strong constraints on mass, leaving 
the posterior close to the prior. At the highest masses 
stochasticity is unimportant, but depending on the as¬ 
sumed metallicity and choice of tracks, the inferred mass 
can vary by up to ^ 0.5 dex. When the distribution is 
dominated by a small number of clusters, this translates 
directly into variations in the overall distribution. Age 
distributions are somewhat less robust that mass ones. 
They are at least mildly dependent on the choice of prior 
at all ages, and at ages < 3 Myr when nebular emission 
contributes significantly to the light, they are very sen¬ 
sitive to metallicity, choice of tracks, and the assumed 
efficiency with which ionizing photons are converted to 
nebular light within the observational aperture. They 
are also quite dependent on the choice of prior at young 
ages. 

The method we develop in this paper is executed in an 
automated pipeline, which is efficient enough that we can 
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derive marginal posterior PDFs for the high confidence 
cluster catalog of 621 clusters in well under an hour on a 
workstation-level machine. This pipeline is available at 
https: //www. slugsps. com, and will form the basis of a 
full stochastic cluster catalog derived from all the star 
clusters in the LEGUS sample. These will be expanded 
to include Ha data from the Ha-LEGUS program (PI: 
R. Ghandar) as they become available. The resulting 
data set will provide an unprecedented sample of cluster 
properties, with realistic uncertainty distributions, from 
across the star-forming part of the Hubble sequence. 
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